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This work seeks to explore and improve the current time-stepping schemes used in 
computational fluid dynamics (CFD) in order to reduce overall computational time. A 
high-order scheme has been developed using a combination of implicit and explicit (IMEX) 
time-stepping Runge-Kutta (RK) schemes which increases numerical stability with respect 
to the time step size, resulting in decreased computational time. The IMEX scheme alone 
does not yield the desired increase in numerical stability, but when used in conjunction 
with an overlapping partitioned (multi-block) domain significant increase in stability is 
observed. To show this, the Overlapping-Partition IMEX (OP IMEX) scheme is applied 
to both one-dimensional (ID) and two-dimensional (2D) problems, the nonlinear viscous 
Burger’s equation and 2D advection equation, respectively. The method uses two differ- 
ent summation by parts (SBP) derivative approximations, second-order and fourth-order 
accurate. The Dirichlet boundary conditions are imposed using the Simultaneous Ap- 
proximation Term (SAT) penalty method. The 6-stage additive Runge-Kutta IMEX time 
integration schemes are fourth-order accurate in time. 

An increase in numerical stability 65 times greater than the fully explicit scheme is 
demonstrated to be achievable with the OP IMEX method applied to ID Burger’s equa- 
tion. Results from the 2D, purely convective, advection equation show stability increases 
on the order of 10 times the explicit scheme using the OP IMEX method. Also, the domain 
partitioning method in this work shows potential for breaking the computational domain 
into manageable sizes such that implicit solutions for full three-dimensional CFD simu- 
lations can be computed using direct solving methods rather than the standard iterative 
methods currently used. 


Nomenclature 


ARK Additive Runge-Kutta 

CFD Computational Fluid Dynamics 

CFL Courant-Friedrichs-Lewy condition 

CDR Convection-Diffusion- Reaction 

DFL Diffusive condition 
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DGFEM 

Discontinuous Galerkin Finite Element Method 

ERIC 

Explicit Runge-Kutta Method 

ESDIRK 

Explicit Singly Diagonally Implicit Runge-Kutta Method 

FD 

Finite-difference 

IMEX 

Implicit-Explicit Method 

IVP 

Initial Value Problem 

OP IMEX 

Overlapping-Partition Implicit-Explicit Method 

ODE 

Ordinary Differential Equation 

NSE 

Navier-Stokes Equations 

PDE 

Partial Differential Equation 

RK 

Runge-Kutta Method 

RK4 

Classic Fourth-order Runge-Kutta Method 

SAT 

Simultaneous Approximation Term 

SBP 

Summation by parts 

ID 

One-dimensional 

2D 

Two-dimensional 

3D 

Three-dimensional 


I. Introduction and Background 

The full three-dimensional (3D) Navier-Stokes equations (NSE), which govern fluid flow, are composed of 
unsteady, convective, and diffusive terms. These equations are nonlinear, due to the convective term, making 
them challenging to solve. Time-dependent computational fluid dynamics (CFD) solutions require numerical 
methods with the ability to efficiently and accurately model each of the mathematical terms, yielding insight 
into the flow physics which cannot or may be very difficult to gain otherwise. 

The general algorithm for numerical integration of time-dependent equations has remained essentially 
the same for the last 3 or 4 decades. 1 Explicit schemes are the most popular and widely used method 
due to their easy implementation and computational efficiency. However, explicit schemes are restricted by 
the Courant-Friedrichs-Lewy (CFL) condition for convective terms and the equivalent diffusive condition 
(DFL) for diffusive terms, which dictate the time step size based on grid size and flow velocity or viscosity, 
respectively. When systems become stiff, the CFL and DFL conditions severely limit the time step size and 
explicit schemes become inefficient. A stiff system, as defined by Chapra and Canale, 2 is a system which 
involves both rapidly and slowly changing components. An example of a geometry-induced stiff system 
important in CFD is that of a boundary layer on an airfoil. The cells near the airfoil wall are very small in 
order to capture the small boundary layer, whereas the cells away from the wall can be orders of magnitude 
larger. One other example resulting in a stiff system is the reaction terms in convection-diffusion-reaction 
(CDR) problems. These examples require solutions to stiff equations which are not efficiently solved using 
explicit schemes. Stiff systems are most efficiently solved using unconditionally stable implicit schemes, 
which are not dependent upon the CFL and DFL conditions. 

One approach to improving the severe time step stability of explicit schemes is to solve individual cells 
or elements using local time-steps, which is often called multi-rate integration. The local time-stepping 
allows for a relaxed stability condition in areas that are less stiff, improving the overall efficiency of the 
method compared to global explicit time-stepping methods. There exist many examples of such schemes, 
some of which will be discussed next. Osher and Sanders 3 introduced a local time stepping method for 
one-dimensional conservation laws. Berger and Oliger 4 present an adaptive mesh refinement method that 
refines the domain in space and time based upon local truncation error. Another adaptive mesh refinement 
algorithm based on cell size and the stability condition is given by Flaherty et al. 5 Dawson and Kirby present 
a first- and second-order accurate time discretization scheme based on local CFL condition, rather than a 
global CFL condition. 6 Tan et al.' present a similar approach using a moving mesh. Bernacki et al. use an 
explicit leap-frog time scheme for use in solving non-dissipative wave propagation problems. 8 While these 
methods improve the explicit stability condition, most of them are implemented at a lower order accuracy, 
second-order accurate in time or less. There exist multi-rate integration methods using high-order accuracy, 
however, they suffer from increased implementation complexity. 

Another approach to improve the time step stability limitations of explicit schemes was originally de- 
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veloped to solve the stiff term in CDR equations implicitly while solving the non-stiff terms explicitly as 
noted by Ascher et al., 9 referred to as Implicit-Explicit (IMEX) methods. They compare various IMEX 
methods and find that they are not a universal cure for all problems. However, they show that IMEX 
schemes can be very effective on certain problems when chosen prudently. Further work from Ascher et al. 10 
focuses strictly on IMEX-RK methods for solving CDR problems, concluding that they exhibit undesirable 
time-step restrictions unless diffusion strongly dominates. Calvo and Novo 11 present two variable step linear 
IMEX schemes for CDR problems, which is competitive with alternative stiffly accurate time integrators. 
Fritzen and Wittekindt 12 show a family of partitioned RK schemes for integration of semi-discrete equations. 
Second- and third-order semi-implicit schemes are used by Zhong 13 for the simulations of hypersonic flows. 
Similar high-order semi-implicit methods are used in Ref. 14. Some limitations of the methods described 
above include: errors in the coupling of the IMEX schemes which lead to less robust simulations, poor 
stability properties, and no error control. 

Kennedy and Carpenter 15 develop additive Runge-Kutta (ARK) methods which allow domain-partitioned 
IMEX time-splitting, meaning that portions of the domain are solved implicitly while other portions are 
solved explicitly. Regions of the domain are implicitly integrated using an L-stable, stifffy-accurate explicit, 
singly diagonally implicit Runge-Kutta (ESDIRK) method, while other parts of the domain are integrated 
using an explicit Runge-Kutta (ERK) method. These are the Implicit-Explicit Runge-Kutta (IMEX-RK) 
schemes that are used in this work. Kanevsky et al. 1 use these additive IMEX-RK schemes to solve systems 
with large amounts of geometry-induced stiffness using a discontinuous Galerkin finite element method 
(DGFEM) discretization. Their work shows that the IMEX-RK method is more efficient than ERK schemes 
when solving geometry-induced stiff systems. They alleviate the time step restriction by solving the very 
small elements in the domain using the ESDIRK scheme while the larger elements are solved using the ERK 
scheme. 

The focus of this work is not to separate the solution of stiff and non-stiff terms as was the focus of 
most the IMEX methods previously described. This work uses the IMEX-RK methods from Kennedy and 
Carpenter 15 to split up the computational domain into partitions or blocks. The independent block partitions 
are advanced in time using the implicit scheme, while the terms coupling the blocks together are advanced 
in time using the explicit scheme. In addition, the partitions are overlapped at the block interfaces, where 
portions of the finite-difference (FD) matrix are solved twice in separate implicit solve blocks. The added 
cost of solving portions of the domain twice are justified by the significant increase in numerical stability, as 
compared to purely explicit schemes, resulting from what will be called the Overlapping-Partition Implicit- 
Explicit (OP IMEX) method. The novel OP IMEX method presented in this work not only yields an increase 
numerical stability, but provides a way by which implicit methods used in CFD can potentially be solved 
using direct solving methods rather than iterative methods. 

II. Methodology 

The numerical methods, theory, and methodology that make up the OP IMEX method will be described 
in this section. This will include discussion of the FD operators formulated using the summation by parts 
(SBP) rule, 16 as well at the Simultaneous Approximation Term (SAT) penalty method of imposing the 
Dirichlet boundary conditions in the weak sense. 17 A general OP IMEX algorithm will also be presented, 
describing how the method was implemented in the code used to produce these results. The partitioning 
and overlapping of the domain is an important aspect of this work and the specific one-dimensional (ID), 
and two-dimensional (2D) domain partitioning will be presented in detail. 

The governing equations of interest in CFD are typically partial differential equations (PDEs), that is, 
equations that are dependent upon multiple variables and their partial derivatives. 18 In order to numerically 
solve such equations, one must discretize them such that the original PDE is transformed into a system 
of algebraic equations. One way to accomplish this is by semi-discretizing the computational domain of 
interest using FD approximations of the spatial derivatives, resulting in a senri-discrete system of ordinary 
differential equations (ODEs). The ODE, by definition, is an equation that is a function of a single variable 
and its derivatives. The semi-discrete scheme is written as an initial value problem (IVP) 

^ = F(f,U(f)), U(i„) = U 0 , (1) 

where U is the unknown variable in vector form of length TV, and N is the number of ODEs resulting from the 
spatial discretization. The system of ODEs are integrated in time using a high-order Runge-Kutta method 
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(see section C). The way in which the boundary conditions are imposed play a vital role when solving PDEs 
and will be discussed in section B. 


A. Finite Difference: Summation by parts 

In this section, the formulation of the FD-based spatial derivative approximations, both first and second 
derivatives, will be developed. The SBP property allows one to obtain a discrete energy estimate, resulting 
in desired stability properties. 


1. First Derivative Approximation 

The domain, D £ [xl,Xr], is discretized using N equidistant points, where x represents the discrete solution 
points, 

X = (xi,x 2 , ..., x n ) t , Xi = x L + (* - l)dx, * = 1,2. ..,7V, dx = X ^ _‘ 7 j L . (2) 

The name summation by parts comes from the fact that the derivative operators mimic integration by parts, 


rXR pxr 

/ <jm x dx = - / (j) x udx. 

J XT, J X Tj 


( 3 ) 


The SBP rule is held by constructing the first derivative approximation, u x (x) = Vu + T P 2 PP , where T P 2 PP 
is the truncation error associated with the order of the scheme, p. The scheme has the following properties: 

V = V~ 1 Q, V = V t 7 ( t V( > 0, C^O, 

(4) 

Q t = B-Q , B = diag(— 1, 0, ..., 0, 1). 


The diagonal V matrix incorporates the grid spacing into the derivative approximation. The Q matrix is 
nearly skew-symmetric, all rows sum to zero, and the first and last columns sum to -1 and 1, respectively. 
One can show the summation by parts property by numerically integrating the following, 

rXR 

/ <j>u x da; « <t> T VVu. (5) 

Jx L 

An inner product of and the approximate derivative Du is formed, where = (</>(a;i), (j)( x 2 ), 4 > { x n)) T 

and the integration weights are in the V norm. Using the properties from equation 4 the summation by 
parts rule is shown as follows, 


<p T VVn = cfVV-'Qu = <t> T (B - Q T ) u 
= 4>nun - 4>iui - 4> t Q t u. 


( 6 ) 


The final form 

<j) T W u = (/>nUn - 4>iUi - <j> T V T Vu, (7) 

<t>U \Z% J** <t>a>udx 

mimics integration by parts shown above in equation 3. 

The truncation error, T P 2 p P , of the approximation has 2 p interior accuracy and p accuracy at the bound- 
ary. The specific SBP approximations used in the work are SBP 1-2-1 and SBP 2-4-2, where the 1-2-1 scheme 
is first-order accurate at the boundary and second-order accurate in the interior, and the 2-4-2 scheme is 
second-order accurate at the boundary and fourth-order accurate in the interior. Refer to appendix A for 
the full SBP operators, which come from Mattsson and Nordstrom. 16 


2. Second Derivative Approximation 


The viscous terms in the equations require a second derivative approximation as well. Integration by parts 
leads to 


<t>u xx dx = <jm x \H 


(j> x u x dx. 


( 8 ) 
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It is possible to use the first derivative approximation to construct the second derivative approximation. 
However, there are several drawbacks to doing this described in, 16 one of which is the width of the operator 
would increase from 2p + 1 to Ap + 1, resulting in a less efficient scheme. The properties of the second 
derivative operator are defined as 

V 2 = V~ 1 {-M + BV), M=M t , C T MC>0, VC, (9) 

where M is restricted to 

M = v T vv + n , n = n, c T ^C>o, vc, (10) 

in order to show the SBP property (see below). Integrating again using the V norm as the integration 
weights yields 

<j) T VV u = <p T VV-H-M + BV ) u = cj) T BV u - cj) T Mu. (11) 

Using the restricted definition of M from equation 10, the final form of the equation is 

(f) T Wu = (fipjUjy — 4>iUi — (j) T V T Wu —c/) T TZu. (12) 

/C“ <t>xU x dx 

The 1Z matrix is known as the remainder matrix which scales with the grid spacing at the same order as 
the error of the second derivative operator. 19 Therefore, the last term of equation 12 is small and decreases 
with increasing resolution so the SBP property is preserved. The SBP second derivative operators used in 
this work can be found in appendix A. 


B. Boundary Conditions: SAT Penalty Method 

The SAT method, coined by Carpenter et ah, 17 is the boundary condition used in conjunction with SBP 
operators. They describe why the SBP property, by itself, is not guaranteed to be stable or time stable. 
The reader is referred to the original paper for an in-depth discussion. In summary, imposing the physical 
boundary condition directly, u{xn) = g(t), can destroy the SBP property and an energy estimate becomes 
unobtainable. Rather, the boundary condition is imposed as an additional ’penalty’ term in the differential 
equation, preserving the SBP property and allowing an energy estimate to be obtained. 

The hyperbolic scalar equation 

= A * e [0,1], t € [0, oo), (13) 

will be used to illustrate the formulation and general theory of the method in this section. For A > 0, the 
wave is traveling to the left and the boundary condition is u(l,t) = g(t). The energy rate is 

u 2 (x,t)dx = A(u 2 (l,f) — u 2 (0,t)), (14) 

and is used with the SBP property to show the necessity of the SAT method for imposing the boundary 
conditions. 

The equation is semi-discretized using the SBP first derivative approximation (see section A) 





(15) 


where u represents the approximate solution at the discrete locations, x, from equation 2. It is important 
to note that the SBP properties from equation 4 give the following energy rate, 


where 


d E 2 2 

= 90, o^o + Qn,nu n , 

(16) 

E(t) = i(u(t),W>u(t)). 

(17) 
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The boundary is imposed as a penalty term as follows, 


= XQu- Tq N}N S(u N - q(t)), 


where S = H 1 (0, 0, ...0, 1) T . Multiplying by the matrix H 

du 


TLV — — XHQu — tqn,n{ 0, 0, ...0, 1 ) t un, 


which can be written as an energy rate 

d E(t) 
d t 


— Q0,0 U 0 + qN,NU% — TqN,NV?N- 


(18) 


(19) 


(20) 


Immediately it is clear that if r > 1 then scheme is both stable and time stable. It can also be proven that 
the SAT method maintains the order of accuracy of SBP approximation (refer to the original paper for the 
proof). 


C. IMEX-RK Time Integration 

The general PDE is reduced to an ODE by discretizing the equation and approximating the spatial derivatives 
using the SBP operators discussed in section A. This results in an ODE at each discrete solution point, and 
therefore, a system of ODEs. Recall that the boundary conditions are included in the differential equation 
using the SAT method (see section B). The semi-discrete equation is now a system of ODEs and is written 
as an initial value problem (I VP) 


— =F(i,U(i)), U(t 0 ) = Uo(ar), (21) 

where U is the solution vector of length N, and N is the number of ODEs corresponding to the number of 
discrete solution points in the domain. The left hand side of the equation, ^ , can now be advanced in time 
by numerical integration, U(t + At) = IT n+1 b There are several forms of numerical integration, however, 
fundamentally they are either explicit, implicit, or some combination of explicit and implicit. This work is 
focused on a newly developed OP IMEX method, which implements high-order ARK schemes by Kennedy 
and Carpenter. 15 

RK schemes can be written in a tabular form, known as the Butcher tableau, 20 as seen below. 


Ci 


a ij 

h 


The Ci vector contains the RK intermediate time levels, the a ij matrix contains the weights for the ith RK 
stage, and the 6; vector contains the integration weights for the final stage. 

The RK schemes used in this work are written in a slightly modified form of the Butcher tableau, 



where bi contains the integration weights for the embedded scheme. The embedded scheme is one order 
less accurate than the main scheme and is used for the adaptive time step controller. The adaptive time 
step controller is used to assess the stability of the 2D OP IMEX method and will be discussed in detail in 
section G. 

This method uses the ARK schemes of Kennedy and Carpenter, 15 which include a 6-stage fourth-order 
accurate ERK scheme used for the explicit time stepping, and a 6-stage fourth-order accurate ESDIRK 
scheme used for implicit time stepping (see appendix A for the schemes in Butcher tableau form). The 
two schemes use the same bi and Cj vectors, such that the intermediate time levels and the final time step 
weights are the same. This is vital for maintaining accuracy when combining the two different schemes for 
the IMEX temporal advancement. Note that the explicit intermediate weights, are zero on and above 
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the main diagonal, af x = 0, j ^ i. The implicit intermediate weights, A]™, are j on the main diagonal, 
A = \,j = i(i > 1), and zeros above the main diagonal, A-™ = 0, j > i. The ESDIRK scheme contains 
explicit contributions to each implicit stage, however, it is important to note that it is an implicit scheme 
and behaves as such. 

Now that the RK schemes used to integrate the system of ODEs in time have been explained, the details 
of how the implicit and explicit contributions to what will be called the ’residual’ will be discussed. The 
residual is referred to the right hand side of equation 21, F(t, U(t)), which is dependent upon explicit and/or 
implicit contributions depending on the location in the computational domain. The intermediate stage values 
are computed using 


i— 1 i— 1 

ub) = u (n) + A a tj F ex (t + CjAt, U&) + A A iT F<m (* + CjAt, U«) 

3 = 1 1=1 ( 22 ) 

+ A tAC F tm {t + Cl At,U {i) ), 

where F ex is the residual with contributions from the explicit terms, F rm is the residual with contributions 
from the implicit terms, is the approximate solution at the intermediate time level, is the approxi- 
mate solution at the previous stage, and At is the chosen time step. The system of ODEs must be advanced 
to each stage simultaneously, meaning that all the implicit and explicit residuals are advanced to the first 
stage before advancing to the second stage and so forth. 

The final stage value is computed using all the stored explicit and implicit residuals for each stage, 

S S 

U (n+1) = U (n) + A t^bi F ex (t + c,At, U (i) ) + AtJ2 h F im (t + aAt, U (i) ). (23) 

i— 1 i= 1 


D. Nonlinear Implicit Solutions 

Newton’s method is implemented in the IMEX method for nonlinear implicit solves. One must calculate the 
residual vector as follows: 


i — 1 

G(U«) = U (i) - U (n) - AtJ2 a tj F ex {t + Cj At, U (j) ) 

3 = 1 

i — 1 

- A/ ^ ,i;;- F im (t + CjAt., U w ) - At Aj™ F im (t, + CiAt, U (i) ); 
3=1 


(24) 


by taking the derivative of equation 24 with respect to the solution stage value, U ( *\ one arrives at the 
Jacobian matrix, 


dG 

<9U« 


M-A tA™ 


QY im 

auib ’ 


(25) 


where [/] is the identity matrix. A system of linear equations is formed with the Jacobian matrix, residual 
vector, and the change in the solution vector as follows, 


-G(U 


(*) \ 
fc+D 


dG 


5 u. 


(26) 


and <5u is obtained by solving the system equations. The change in the solution vector, Ju, is used to update 
the velocity vector, 


U 


M 

fc+i 



+ <5u, 


(27) 


and Newton’s method is iterated until convergence is reached. The number of Newton iterations required to 
obtain convergence is dependent upon the nonlinearity of the system of equations and the desired convergence 
tolerance. This solution procedure is repeated for each stage in the RK scheme. 

As a side note, this nonlinear system of equations solution procedure can also be applied to linear 
systems of equations, where Newton’s method will converge in one iteration for the linear set of equations. 
In this work, the nonlinear Burger’s equation and linear 2D advection equation were studied using the same 
nonlinear implicit solution procedure described above. 
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E. ID Domain Partitioning and Overlapping 

The domain partitioning and overlapping used in the OP IMEX method, applied to a ID problem, is 
presented in this section. This will include details of the method using both SBP 1-2-1 and SBP 2-4-2 
schemes, respectively. In this work, the computational domain x £ \xl,xr\ is split into three evenly spaced 
partitions: xi £ [xl,xa], x 2 £ [xa+i,xb], x 3 £ [xb+i,xr\. The number of solution points is adjustable, 
however, it must be divisible by three in order to always maintain the three equispaced partitions. The 
choice of three equispaced partitions is not required by the OP IMEX methodology and was chosen simply 
for easier code implementation for this initial study. 


1 . SBP 1 - 2-1 Scheme 

The second-order spatially accurate scheme (SBP 1-2-1) uses a three-point wide stencil, requiring the so- 
lution of a tridiagonal matrix for implicit time-stepping. Splitting the FD matrix into partitions or blocks 
results in portions of the domain that can be independently solved, which is desirable when running large 
simulations on multiple processors. The splitting and overlapping of the FD matrix is shown with the fol- 
lowing simple example. The FD matrix results from a computational domain with eight solution points, 
a = (01,02, ...,as) T , which are along the main diagonal, 


a 1 ci 

&2 <22 C2 

^3 <23 C 3 

64 CI4 

C 4 

b 3 

a 5 c 5 
&6 O'd ^6 

67 dj Cj 

b$ as 


(28) 


The domain is partitioned at a single location, between 04 and 05 as shown. The upper left and lower right 
blocks are implicitly solved (shown inside red blocks) while C4 and 65 (outside red blocks) are explicitly 
solved. The simple partitioning of the FD matrix above, along with the implicit and explicit solves, does not 
gain much in stability by itself. The real gain in stability comes from overlapping the partitions and solving 
small portions of the domain twice. The system of equations, shown in equation 28, will be used to show the 
structure of the FD matrix as the partitions are overlapped. The computational domain a = (ai, 02, ..., a 3 ) T 
is again partitioned between solution points 04 and 05. However, this time the partition is overlapped once 
as follows, 


ai ci 

b 2 <22 c 2 
b 3 a 3 c 3 
64 04 C4 
65 <25 

C5 

64 

a 4 C4 
^5 C5 

&6 Cq 

bj dY Cy 
bg ag 


(29) 


The system of equations is now 10x10 due to the single overlap of each partition, rather than 8x8. The 
solution points <24 and <25 are solved in the upper left block and also in the lower right block. The explicit 
contributions are shifted to C5 and 64 which are also included in the implicit block solves as seen above. 
This shifted matrix structure, resulting from the overlapped partitions, is mathematically demonstrated 
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using a bordering technique and matrix operations or visually by reordering the domain and applying FD 
discretization. The partitions can be overlapped multiple times; the explicit contributions continue to shift 
outward. In this final example of the tridiagonal FD matrix with eight solutions points, the matrix is 
partitioned again between 04 and a 5 and overlapped twice, 


ai ci 
^2 0-2 C2 
&3 «3 C3 

64 a 4 C4 

65 05 C5 
be ae 

ce 

be 

CL3 C 3 
64 CL4 C4 

65 a 3 C5 

&6 Qd 

&7 dj C7 
&8 &8 


(30) 


resulting in a 12x12 system of equations. The solution points 03, 04, as, and a6 are implicitly solved (and 
their neighboring contributions) in both the upper left and lower right blocks (shown in red boxes). The 
explicit contributions have moved out and are now C6 and 63. Overlapping the partitions increases the 
computational expense of the solution by solving larger set of equations, however, if the gain in stability 
outweighs the added cost then the method is more efficient. The simple example used above would be very 
inefficient because a majority of the solution points are double solved with just two overlaps, nearly doubling 
the size of the original system of equations. However, when solving large systems of equations the overlapped 
duplicate solves are only a small fraction of the total number of solution points and the method has the 
potential of being very efficient. 


2. SBP 2-^-2 Scheme 

The fourth-order spatially accurate scheme (SBP 2-4-2) has a five-point wide stencil and results in a penta- 
diagonal FD matrix. The domain partitioning and overlapping is performed in a similar manner as with the 
SBP 1-2-1 scheme. However, the added diagonals in the matrix change the explicit contributions. To show 
this, a FD matrix with ten solution points, a = (ai,a2, ...,aio) T , is shown 



(31) 


The domain is split between solution points 05 and a e and the implicit solves are shown with the red boxes. 
The wider stencil results in more ghost points that are explicitly solved: e4, C5, and e$ for the upper left 
block, and de, be, and dj for the lower right block. A similar behavior in the shifted matrix, as compared 
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to the narrow SBP 1-2-1 scheme, is seen when overlapping the partitions. A single overlap results in this 
matrix structure 



\ 


(32) 


Again, the implicit solves are shown inside the red boxes and the explicit contributions (e$, C 6 , e 6 for the 
upper left block and 65 , de for the lower right block) have shifted outward. Note that with a single 
overlap, duplicate explicit contributions are found inside the implicit solve blocks for both the tridiagonal 
and pentadiagonal matrices. A final example of a pentadiagonal matrix with two overlaps, 



(33) 


As expected, the explicit contributions have shifted two more points outward. It is important to remark 
that the SBP operator matrices do not form perfect tridiagonal and pentadiagonal matrices as shown in the 
above examples. The upper left and lower right portions of the matrices are modified to handle boundary 
accuracy (see appendix A). The matrix partitions must be sufficiently away from the boundaries so as to 
avoid the boundary modified portion of the matrix. 


F. 2D Domain Partitioning and Overlapping 

The 2D domain partitioning and overlapping, for the OP IMEX method, is presented in this section. A 
simple example is presented, showing how the domain is partitioned and mapped to the full matrix for the 
implicit temporal advancement. 

The following example demonstrates the domain ordering and mapping of a 2D domain discretized into 
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a 6 x 6 matrix, 


25 

26 

27 

34 

35 

36 

22 

23 

24 

31 

32 

33 

19 

20 

21 

28 

29 

30 

07 

08 

09 

16 

17 

18 

04 

05 

06 

13 

14 

15 

01 

02 

03 

10 

11 

12 


( 34 ) 


The domain is split into four blocks, each containing a 3x3 portion of the domain. The block ordering 
goes from left-to-right, bottom-to-top, as does the individual block ordering of the solution points. In the 
example above (equation 34): block 1 is outlined in blue and contains grid points 1-9, block 2 is outlined in 
yellow and contains grid points 10-18, block 3 is outlined in red and contains grid points 19-27, and block 4 
is outlined in green and contains grid points 28-36. The ordering of the domain is such that the grid points 
in each block are clustered together, which yields smaller implicit block solves once mapped to the full FD 
matrix (see below). 

In order to implicitly solve the 2D problem, a system of equations is formed containing an equation for 
each grid point in the domain. The example domain (6x6) from above is mapped to a full 36x36 matrix, 
which is shown next. 



19 20 22 

19 20 21 23 

20 21 24 

19 22 23 25 

20 22 23 24 26 

21 23 24 27 

22 25 26 

23 25 26 27 

24 26 27 


21 


34 


28 29 31 

28 29 30 32 

29 30 33 

28 31 32 34 

29 31 32 33 35 

30 32 33 36 

31 34 35 

32 34 35 36 

33 35 36 


(35) 


The colored boxes in equation 35 correspond to the block number from equation 34 and represent the 
implicitly solved partitions (block 1 is blue, block 2 is yellow, block 3 is red, and block 4 is green). The terms 
outside the boxes are the overlapping ghost points which are advanced explicitly. Note that the re-ordering 
of each block from above yields implicit solve blocks that are tightly clustered together. 

G. Time Step Controller 

Evaluating the stability of the 2D advection equation was less cut-and-dry as compared to the ID nonlinear 
Burger’s equation. The stability envelop of the ID simulation of Burger’s equation was analyzed simply by 
adjusting the time step and running the code. As the time step became too large to maintain numerical 
stability, the simulation would diverge very quickly due to the nonlinear nature of the equation. The 2D 
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advection equation, on the other hand, did not demonstrate a clear cut stability envelop. As the time step 
was increased, the error in the solution increased but it did not quickly diverge, making it difficult to quantify 
the stability of the method. This was most likely due to the linear behavior of the equation. Therefore, 
an adaptive time step controller was implemented to ’automatically’ quantify the stability of the OP IMEX 
method on this 2D problem. 

The explanation from Kanevsky et al. 1 is reproduced here, giving a nice summary on time step controllers. 
The embedded time-integration scheme, which is one order lower than the main scheme, is used to compute 
a temporal error. The temporal error is then passed into the adaptive time step controller, such as an I, PI, 
or PID controller, which automatically and adaptively adjusts the size of the time step. 

The I-based controller derivation is shown next, giving insight into the general theory of adaptive time 
step controllers. The temporal error, S , is computed by subtracting the solution based on the embedded 
scheme of order p from the main scheme of order p + 1 

<5 = U — U 

= (U exact + 0((At) p+1 )) - (U 

exact + 0((At) p )) _ 

= 0({At) p ) 

= C(At) p , 


where C is some constant. The temporal error between two different time steps <5^ and is compared 


$("+!) _ C((At)( n+1 ))P 
<5( n ) - C((A t) n )p 
/(A t)("+ 1 >\ p 

~V (At)" ) 


(37) 


where (At)(" +1 ) is the time step to be computed based on a temporal error of The temporal error is 

then given as a user-defined input, e = m n+1 ), and equation 37 is rearranged to solve for the computed time 
step, 


(Af) (n+1) 



(38) 


Lastly, a factor of safety k is added and the computed temporal error is written as an L m 
controller is 


(At)| n+1) 


k(A t) n 



norm. The I-based 
(39) 


The PID-based controller uses the computed temporal error from two previous time steps, in addition to the 
current time step, and goes as follows, 


(At) 


(n+l) 

PID 



(ll^- 2) ll 


A 


(40) 


where e is the specified tolerance for the temporal error, and p is the order of accuracy of the embedded 
scheme. The temporal error is computed as 


where 


£(«+!) = u( n+1 ) - 

S S 

Uffi+i) = A t^bi F“’(U (,; )) + A tJ2 b < F im (u (i) ), 

i=l i= 1 

s s 

uffi+i) = At^6,F ea ’(U (i;) ) + Af^6 2 ;F im (U (i) ). 
Simplifying the above equations, 


i = 1 


i=l 


<y( n+1 > = At ^2(bi - Si) F ea: (U (i) ) + A t^2(bi - hi) F™(U (i )). 

i=l i= 1 


(41) 


(42) 


(43) 
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The PID controller for this work uses the gains taken from Kennedy and Carpenter, 15 


kj = 0.25, k P = 0.14, kp = 0.10, u) n = 1, 


(44) 


where 


Therefore, 


pa = [ki + kp + 



pfi = [kp + 2 uj n k D ], 


pX 


2 ui 


kp. 


.49 a .34 , .10 

a = — , p = — , A = — . 
P P P 


(44) 

(44) 


Also, an L 2 norm of the computed temporal error is used rather than an L ^ norm of the temporal error. 
Refer to 23-25 for more information on adaptive time step controllers. 


H. Overlapping-Partition IMEX Algorithm 

The OP IMEX-RK algorithm is presented in this section to further clarify the implementation of the method. 
The discrete solution vector, u, is passed into the routine and advanced in time by one time step as shown 
here. The routine is sequentially repeated to advance the solution to any given final time. The algorithm 
loops through all the block partitions, advancing the residuals containing the implicit contributions (F* m ) 
to the first RK stage. After which, the algorithm loops through all the block partitions again, advancing 
the residuals containing the explicit contributions (F ea: ) to the first RK stage. This is repeated and the 
discrete solution vector is advanced through all six RK stages such that the implicit and explicit residuals 
are advanced simultaneously. The linear solver used for the implicit solutions of the ID problem is a freely 
distributed software called SPARSKIT. 22 Due to the sparse matrices resulting from the method, a direct 
linear solver tailored for sparse matrices was used (iluk and lusol routines specifically). The 2D problem 
required additional mapping of the computational domain to the solution vector and required a different 
sparse matrix solver. PETSc library 26 was used for the implicit linear solutions, specifically the PCLU 
routine from the KSP solvers, and gave good results. See algorithm 1 below for a pseudo code of the scheme. 

Algorithm 1 Overlapping-Partition Implicit-Explicit Runge-Kutta Algorithm 
for i = 1 — > s do ! loop over s-stages 

for e = 1 — >• nparts do ! loop through all partitions 

for j =— ; > i — 1 do ! explicit integration of implicit scheme 
u = u n + At * afj * F ex + A t* Alj 1 * F res 
end for 

iihat — u ! store total explicit residual 

while ( newton res < 1.0 E — 012) do ! implicit solve (Newton’s method) 

Call implicit residual/ Jacobian routine (solve for F lm and (^/) 

G res = Uhat — u + At * A lm * F lm ! Newton’s residual 
^Su= —G res ! Newton’s method 
newton res = \J ( G res ) 2 
Call linear solver (solve for 6u) 
u = u + Su ! update u 
end while 
end for 

for e = 1 — > nparts do ! loop through all partitions 

Call explicit residual routine (solve for F ex ) 

end for 
end for 

! final integration using RK weights (b(k) vector) 

for k = 1 — > s do 

u = u n + At * b(k) * F ex + A t* b(k) * F irn 

end for 
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III. ID Numerical Test: Burger’s Equation 

The nonlinear viscous Burger’s equation written in conservative form is 

= tu xx , (44) 

where e is the kinematic viscosity, and u is the velocity of the fluid. Burger’s equation is a model equation 
of the NSE due to its unsteady term, nonlinear convective term, and diffusive term, which are much like the 
terms found in the NSE. Applying new numerical methods to Burger’s equation is a common way to gain 
understanding of how a method performs on nonlinear hyperbolic PDEs, without the added complexities of 
two- and three-dimensional equations. The newly developed OP IMEX method, presented in section II, is 
first applied to and tested on Burger’s equation. Details of the implementation and the ID OP IMEX results 
will follow in this section. 


A. Semi-discrete Equation and Boundary Conditions 


In order to numerically solve Burger’s equation, the convective and diffusive terms are discretized and ap- 
proximated over the ID computational domain using SBP operators (see section A). The Dirichlet boundary 
conditions are imposed using the SAT penalty method (see section B). 

The semi-discrete Burger’s equation, including SAT boundary conditions and initial condition (u 0 ), is 
written as follows: 

u t = --VUu + eV 2 u + a 0 V~ 1 e 0 [a 0 'Ui - e(<Su)i - g 0 {t)\ 

-I- a\P~ 1 ei [aiUN — e(<Su )n , 


where 

x G [x L ,x R ], 

v(x L ,t) + \v(x L ,t)\ 

go{t) = g v(x L ,t) - ev x {x L ,t), 

v(x R ,t) - \v(x R ,t)\ 

9i{t) = g v{x R ,t) - ev x (x R ,t), 

U = diag(u), e 0 = (1,0, ...) T , e x = (..., 0, 1) T , 
Ul + |ui| _ UN — |ttjv| 
a ° — 3 ’ ai “ 3 ’ 

u 0 = v(x, 0), 


and S is the boundary first derivative operator, 7? is the SBP first derivative operator matrix, D 2 is the 
SBP second derivative operator matrix, and v is the exact analytical solution. Note that the convective 
SAT penalty contributions are formulated to act as a switch: waves convecting into the domain turn the 
switch ’on’ and the penalty is imposed, and waves convecting out of the domain turn the switch ’off’ and the 
penalty is zero. This is observed in do and 01 , and in the first terms of go(t) and g\{t). The energy analysis 
is used to determine that <To = —1 and ay = 1, which is derived by Fisher et al. 21 and will not be presented 
here. The semi-discrete equation is advanced in time at discrete locations, t G [0,T], using the IMEX-RK 
time integration schemes. 


B. Test Problem 

The exact solution to Burger’s equation used to test the ID OP IMEX-RK method is 

v(x, t) = 1 — tanh , 2 ; G [— 10, 10], t € [0,T], (43) 

where e is the kinematic viscosity, Xq = 0, and the initial condition is uq = v(x, 0). The problem is integrated 
up to time, T = 5. The wave propagates left-to-right as the solution advances in time (see figure 1). 
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Figure 1. Discrete solution to Burger’s equation on a grid with 513 equispaced points, using fourth-order 
accurate OP IMEX-RK integration, and fourth-order spatially accurate SBP 2-4-2 scheme. 


1. Grid, Resolution Study 

The following tables show the effects of overlapping the partitions on coarse, medium, and fine grids. The 
fully explicit (no partitions) schemes were run with the maximum time step based on numerical stability, 
which was determined visually, and their corresponding maximum CFL and DFL are presented. Similarly, 
the OP IMEX schemes were run with varying number of overlapped partitions. The maximum CFL and 
DFL of the OP IMEX schemes, based on stability, are presented. Results are shown for both the second- 
order accurate operators (SBP 1-2-1) and the fourth-order accurate operators (SBP 2-4-2) so as to gain an 
understanding of the method’s performance using FD stencils of different widths. The factor of increase 
quantifies the amount of increase in CFL and DFL that is gained by the OP IMEX scheme with respect to 
the fully explicit scheme. 


Table 1. OP IMEX stability results on coarse grid (258 points) compared to the explicit scheme. 


N = 258 



SBP 1-2-1 

SBP 2-4-2 


CFL 

DFL 

Factor of Increase 

CFL 

DFL 

Factor of Increase 

Fully Explicit 

0.6 

1.0 

1.0 

0.5 

0.7 

1.0 

IMEX, overlap = 0 

3.6 

5.8 

5.6 

2.5 

4.0 

5.4 

IMEX, overlap = 1 

12.9 

20.6 

20.0 

10.2 

16.4 

22.1 

IMEX, overlap = 2 

13.8 

22.1 

21.4 

10.3 

16.5 

22.2 

IMEX, overlap = 5 

13.9 

22.3 

21.6 

10.3 

16.5 

22.2 

IMEX, overlap = 10 

13.9 

22.3 

21.6 

10.3 

16.5 

22.2 


N = 513 



SBP 1-2-1 

SBP 2-4-2 


CFL 

DFL 

Factor of Increase 

CFL 

DFL 

Factor of Increase 

Fully Explicit 

0.3 

1.0 

1.0 

0.2 

0.7 

1.0 

IMEX, overlap = 0 

1.5 

4.8 

4.8 

1.1 

3.4 

5.3 

IMEX, overlap = 1 

12.4 

39.7 

40.3 

9.2 

29.5 

45.0 

IMEX, overlap = 2 

12.8 

41.0 

41.7 

9.2 

29.5 

45.0 

IMEX, overlap = 5 

12.8 

40.8 

41.5 

9.2 

29.5 

45.0 

IMEX, overlap = 10 

12.8 

41.0 

41.7 

9.2 

29.5 

45.0 


N = 1014 



SBP 1-2-1 

SBP 2-4-2 


CFL 

DFL 

Factor of Increase 

CFL 

DFL 

Factor of Increase 

Fully Explicit 

0.2 

1.0 

1.0 

0.1 

0.8 

1.0 

IMEX, overlap = 0 

0.7 

4.5 

4.4 

0.5 

3.3 

4.3 

IMEX, overlap = 1 

8.4 

53.2 

51.9 

7.3 

46.2 

60.0 

IMEX, overlap = 2 

10.6 

67.3 

65.6 

7.9 

50.0 

65.0 

IMEX, overlap = 5 

10.6 

67.3 

65.6 

7.9 

50.0 

65.0 

IMEX, overlap = 10 

10.6 

67.3 

65.6 

7.9 

50.0 

65.0 
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2. Viscosity Study 


A study was performed on the affect of overlapping the partitions while varying the viscosity. Similar to 
the grid resolution study, the maximum CFL and DFL are presented for the fully explicit scheme and the 
various overlaps with the IMEX scheme. 


Table 2. OP IMEX stability results on medium grid (513 points) with a diffusion coefficient of 0.0625 


e = 0.0625 



SBP 1-2-1 

SBP 2-4-2 


CFL 

DFL 

Factor of Increase 

CFL 

DFL 

Factor of Increase 

Fully Explicit 

1.2 

0.9 

1.0 

1.0 

0.7 

1.0 

IMEX, overlap = 0 

9.9 

7.9 

8.2 

6.7 

5.4 

6.7 

IMEX, overlap = 1 

14.0 

11.2 

11.7 

11.2 

9.0 

11.2 

IMEX, overlap = 2 

13.9 

11.1 

11.6 

11.1 

8.9 

11.1 

IMEX, overlap = 5 

13.8 

11.0 

11.5 

10.9 

8.7 

10.9 

IMEX, overlap = 10 

14.1 

11.3 

11.7 

11.0 

8.8 

11.0 


e = 0.125 



SBP 1-2-1 

SBP 2-4-2 


CFL 

DFL 

Factor of Increase 

CFL 

DFL 

Factor of Increase 

Fully Explicit 

0.6 

1.0 

1.0 

0.5 

0.7 

1.0 

IMEX, overlap = 0 

3.5 

5.6 

5.7 

0.5 

0.7 

1.0 

IMEX, overlap = 1 

13.2 

21.1 

21.5 

10.2 

16.3 

22.1 

IMEX, overlap = 2 

13.8 

22.0 

22.4 

10.2 

16.4 

22.2 

IMEX, overlap = 5 

13.9 

22.2 

22.6 

10.2 

16.4 

22.2 

IMEX, overlap = 10 

13.9 

22.3 

22.7 

10.2 

16.4 

22.2 


e = 0.25 



SBP 1-2-1 

SBP 2-4-2 


CFL 

DFL 

Factor of Increase 

CFL 

DFL 

Factor of Increase 

Fully Explicit 

0.3 

1.0 

1.0 

0.2 

0.7 

1.0 

IMEX, overlap = 0 

1.5 

4.8 

4.8 

1.1 

3.4 

5.3 

IMEX, overlap = 1 

12.4 

39.7 

40.3 

9.2 

29.5 

45.0 

IMEX, overlap = 2 

12.8 

40.0 

41.7 

9.2 

29.5 

45.0 

IMEX, overlap = 5 

12.8 

40.8 

41.5 

9.2 

29.5 

45.0 

IMEX, overlap = 10 

12.8 

41.0 

41.7 

9.2 

29.5 

45.0 


C. Discussion 

It is observed from both studies that the DFL is the limiting condition (DFL «1 for the explicit schemes 
while the CFL varies) in this problem. The OP IMEX scheme shows promising results; a maximum factor of 
increase of «65 on the fine grid, and «22 on the coarse grid. As the grid is refined (Table 1) the numerical 
stability increases. The viscosity study (Table 2) shows that added viscosity also increases the stability. 
Also, a single overlap in the partitions yields essentially the same increase in stability as multiple overlaps. 
Referring back to the formulation of the overlapping partitions, section E, a single overlap causes the explicit 
contributions to be duplicates of terms contained in the implicit blocks. This appears to be key in the added 
stability of the OP IMEX scheme. 

D. Design Order: Burger’s equation 

Verifying the design order is an important validation of the mathematical correctness of the OP IMEX 
method. The spatial and temporal convergence tests are presented in Figure 2, showing that the OP IMEX 
methodology and the implementation are correct for Burger’s equation. The design order for the OP IMEX 
method with no overlaps and with a single overlap, verifying that overlapping the partitions in this way is 
also correct. All plots are in log-log scale in order to more readily see the order of convergence (as a slope 
rather than a power). Tables containing the exact design order error values are given in appendix A. 

IV. 2D Numerical Test: Advection Equation 

Significant increase in numerical stability was demonstrated using the OP IMEX method applied to a 
ID problem (Burger’s equation). The FD matrix structure for solutions to 2D problems is different than 
those of ID problems. Hence, the next part of this work is to develop the OP IMEX method applied to 2D 
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problems, and to test its performance on a simple problem. The 2D advection equation was chosen for this 
study; it is a linear first-order hyperbolic PDE, 

Ut + A u x + B u y = 0, (43) 

where A is the wave speed in the x-direction, and B is the wave speed in the y-direction. Note that it lacks a 
diffusive term and is purely convective, meaning that only the CFL condition will be studied in this problem 
(no DFL condition). 


A. Semi-discrete Equation and Boundary Conditions 


In order to numerically solve the 2D advection equation, the differential equation is transformed into a system 
of algebraic equations by the following: the domain is discretized, the derivatives are approximated using 
SBP operators, and the boundary conditions are imposed using the SAT penalty method. The derivative 
approximations are computed using SBP operators in a line-by-line fashion (i.e., by looping over the grid 
points in y and performing the derivatives with respect to x and then looping over the grid points in 
x performing the derivatives with respect to y). This is necessary because the SBP operator matrices 
approximate the derivative of a vector or line and is not applicable to a full 2D domain by a single matrix 
operation. 

The senri-discrete 2D advection equation, including SAT boundary conditions and initial condition (uo), 
is written as follows: 

u t = — AVu(:,yi) - BVu(xi , :) 

+ cr 0 'P~ 1 e Q [a L u(xi,yi) - g L (t)\ 

+ a{P~ l ei [a R u(x N , yi) - g B {t)\ (43) 

+ a 0 'P~ 1 e 0 [d B u(xi, yi) - g B (t)] 

+ <Ji'P~ 1 ei [dr u(xi, y N ) - grit )} , 


where 


xe[x L ,x R ], y e[y B ,yT], 

v(x 1 ,y i ,t) + \v(x 1 ,y i ,t)\ 
gUt) = v{x 1 ,y i ,t), 

v(x N ,yi,t) - \v(x N ,yi,t)\ , ,, 

gR{t) = v(x N ,yi,t), 

v(x i ,y 1 ,t) + \v(x i ,y 1 ,t)\ 
gB(t) = v{Xi,yx,t), 

v(xi,y N ,t)-\v(xi,y N ,t)\ 
grit) = v{Xi,y N ,t), 


e o = (1, 0, ...) T , ei = (..., 0,1) 5 


aL = 


A + |A| 


a B = 


A-\A\ 


B+\B\ 


a B = 


a _ r 


u 0 = «(x,y,0), 


B~\B\ 

2 


and 1 > * > N, T> is the SBP first derivative operator matrix, <7o = — 1 and a i = 1 (determined from energy 
analysis), and v is the exact analytical solution. The SAT penalty boundary conditions are the same as 
those imposed for Burger’s equation without the viscous terms. 


B. Test Problem 

The exact solution to the 2D advection equation used to test the 2D OP IMEX-RK method is 

v(x, y, t) = sin(27r(a; — Af))sin(27r(y — Bt)), x £ [0,0.5], y£ [0,0.5], f S [0,T], (42) 

where A is the wave speed in the x-direction, B is the wave speed in the y-direction, and the initial condition 
is given by uq = v(x, y,0). The problem is integrated up to time, T = 7. In figure 3 below, the propagation 
of the 2D wave is shown up to time, T « 2. The wave repeatedly propagates through the domain as the 
equation is integrated farther in time. Note that the magnitude is scaled by 0.1 for better visualization. 
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C. Results 


The numerical stability of the OP IMEX method solving the 2D advection equation was assessed on four 
different grid resolutions: 30x30, 60x60, 90x90, and 120x120. These solution domains are mapped to 
matrices of sizes: 400x400, 3600 x 3600, 8100x8100, and 14400x14400, for the implicit solves, respectively. 
Recall that the stability of Burger’s equation was dependent upon the grid resolution, hence the 2D results 
are given on several grid sizes. Note that there is no diffusive term in the 2D advection equation as it is 
purely convective, meaning that there is no DFL stability criteria. The CFL stability will be discussed in 
this section. 

The variable time stepping error tolerance was varied to determine the stability of the OP IMEX method. 
The linear behavior of the advection equation yielded slowly diverging solutions, unlike the quickly diverging 
solutions of the nonlinear Burger’s equation. This made is necessary to choose a stability parameter, which 
was chosen to be the magnitude of the advecting wave, C, where the exact solution wave magnitude is C = 1. 
The approximate solution was considered unstable when C > 1.5. The edge of the explicit stability regime 
was found by adjusting the error tolerance to be as large as possible while still maintaining a stable explicit 
solution. This was done for both SBP 1-2-1 and SBP 2-4-2 schemes on all four grid resolutions and is plotted 
in figure 4. The CFL oscillates until the controller locks into the stable time step after approximately 100 
and 150 time steps for the SBP 1-2-1 and SBP 2-4-2 schemes, respectively. The second-order scheme (SBP 
1-2-1) is stable at CFL ss 2 and the fourth-order scheme (SBP 2-4-2) is stable at CFL ss 1.5, which is the 
baseline explicit stability that the OP IMEX scheme will be compared against. The L 2 norms are computed 
for the explicit solutions for comparison with the OP IMEX solutions and will be presented later. 

The partitioning of the computational domain into blocks introduces more error into the solution as 
compared to the full FD domain used for the purely explicit baseline runs. A very small error tolerance in 
the controller, on the order of IE-6 or smaller, will drive the OP IMEX solutions to have a smaller CFL 
than the explicit baseline result because of this added error, even though the scheme is numerically stable at 
much higher CFL values. In order to use the controller as a stability evaluation tool, two different controller 
tolerance values, e = 5E-3 and e = IE-2, were used to test the scheme on the different grid sizes. Testing 
showed that those values gave good indication of the upper stability limit of the OP IMEX method. 

Figures 5-8 show how the controller adjusted the CFL over the span of the integration for the explicit 
(EX), implicit (IM), IMEX with no overlaps (IMEX-0), IMEX with one overlap (IMEX-1), IMEX with 
five overlaps (IMEX-5), and IMEX with nine overlaps (IMEX-9) schemes while using an error tolerance of 
e = 5E-3. The CFL of the implicit scheme continues to grow through out the entire integration as it is 
unconditionally stable. The CFL of explicit solution stays at the baseline values of « 2 and ss 1.5 for grids 
30x30 and 60x60, and is unstable for the 90x90 and 120x120 grids. The CFL of the IMEX scheme with 
no overlaps and one overlap is unstable for all the grids expect for the 30x30, where it increases and then 
slowly decreases due to instability. The CFL of the IMEX scheme with nine overlaps is stable for all grids 
and reached a maximum of CFL ss 8 for the SBP 1-2-1 scheme, and CFL ss 7 for the SBP 2-4-2, both on 
the 120x120 grid. A general trend is seen of a slight increase in the CFL as the grid resolution is increased 
for both the 1-2-1 scheme and the wider 2-4-2 scheme. 

Figures 9-12 show the CFL vs. Time Step results with an increased controller error tolerance of e = 
IE-2. The results are similar to those presented in figures 5 - 8 for the error tolerance of e = 5E-3, except 
the increased error tolerance allows the controller to show higher stable CFL values for the IMEX schemes; 
CFL « 10 for the 1-2-1 scheme and CFL ss 9 for the 2-4-2 scheme, when solved on the fine grid (120x120). 
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N 


N 


(a) Spatial design order of two: No overlaps (b) Spatial design order of two: Single over- 
lap 




N N 


(c) Spatial design order of four: No overlaps (d) Spatial design order of four: Single over- 
lap 




(e) Temporal design order of four: No over- (f) Temporal design order of four: Single 
laps overlap 


Figure 2. Burger’s equation L 2 and L 0 0 spatial and temporal design order for OP IMEX scheme. 
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EX 242 _60_0 


EX 242 _60_0 


EX 242 60 0 





(a) t=0 


(b) t=0.21 


(c) t=0.51 





(d) t=0.75 


(e) t=0.99 


(f) t=1.26 





(g) t=1.5 


(h) t=1.74 


(i) t=2.01 


Figure 3. The approximate solution to the 2D advection equation using a 60 x 60 domain, explicit time stepping, 
and SBP 2-4-2 derivative approximations. 
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(a) SBP 1-2-1 Scheme 


(b) SBP 2-4-2 Scheme 


Figure 4. Behavior of CFL vs. Time Step as the adaptive time step controller locks into a stable CFL for the 
ERK schemes on various grid resolutions. 




(a) SBP 1-2-1 Scheme 


(b) SBP 2-4-2 Scheme 


Figure 5. CFL vs. Time step with an error tolerance of e = 5E-3 on the 30x30 grid. 
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CFL CFL 




(a) SBP 1-2-1 Scheme 


(b) SBP 2-4-2 Scheme 


Figure 6. CFL vs. Time step with an error tolerance of e = 5E-3 on the 60x60 grid. 




- EX 


IM 

-V 

IMEX 0 

f — 

- IMEX 1 

<► — 

IMEX 5 


IMEX 9 





20 30 40 

Time Step 


(a) SBP 1-2-1 Scheme 


(b) SBP 2-4-2 Scheme 


Figure 7. CFL vs. Time step with an error tolerance of e 


5E-3 on the 90x90 grid. 
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CFL CFL 




(a) SBP 1-2-1 Scheme 


(b) SBP 2-4-2 Scheme 


Figure 8. CFL vs. Time step with an error tolerance of e = 5E-3 on the 120x120 grid. 




(a) SBP 1-2-1 Scheme (b) SBP 2-4-2 Scheme 

Figure 9. CFL vs. Time step with an error tolerance of e = IE-2 on the 30x30 grid. 
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CFL CFL 




(a) SBP 1-2-1 Scheme (b) SBP 2-4-2 Scheme 

Figure 10. CFL vs. Time step with an error tolerance of e = IE-2 on the 60x60 grid. 




(a) SBP 1-2-1 Scheme 


(b) SBP 2-4-2 Scheme 


Figure 11. 


CFL vs. Time step with an error tolerance of e = IE-2 on the 90x90 grid. 
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CFL 




(a) SBP 1-2-1 Scheme (b) SBP 2-4-2 Scheme 

Figure 12. CFL vs. Time step with an error tolerance of e = IE-2 on the 120x120 grid. 
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The final set of results presented here are plots of the CFL vs. Time Step of the IMEX schemes ran with 
a larger time step error tolerance, e = 5E-2. At this larger error tolerance, the controller allows the time step 
to increase to the upper bound of stability. In fact, several overlaps are required to obtain a stable solution 
with this large error tolerance in some cases. 

Figure 13 shows results of the two schemes on the 30x30 grid. The 2-4-2 scheme was not stable with 1 
or 2 overlaps and reached a CFL« 7 with 4 overlaps. The 1-2-1 scheme was stable at a slightly higher CFL, 
around 8, and was stable with a single overlap. 




(a) SBP 1-2-1 Scheme (b) SBP 2-4-2 Scheme 

Figure 13. CFL vs. Time step with an error tolerance of e = 5E-2 on the 30x30 grid. 


The stability on the 60x60 grid, figure 14, increased from the 30x30 grid but required more overlaps to 
maintain stability; 4 for the 1-2-1 scheme and 6 for the 2-4-2 scheme. The 1-2-1 scheme with 7 overlaps went 
unstable just before the end and is not shown in the plot. 




(a) SBP 1-2-1 Scheme 


(b) SBP 2-4-2 Scheme 


Figure 14. CFL vs. Time step with an error tolerance of f = 5E-2 on the 60x60 grid. 


In figure 15, the 90x90 grid required 7 overlaps for the 1-2-1 scheme, and 9 overlaps for the 2-4-2 scheme 
in order to maintain stability. The results show the same trends: an increase in stability for both schemes as 
compared to the coarser grids, and the 1-2-1 scheme has a slightly higher stability than the 2-4-2 scheme. 
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(a) SBP 1-2-1 Scheme 


(b) SBP 2-4-2 Scheme 


Figure 15. CFL vs. Time step with an error tolerance of f = 5E-2 on the 90x90 grid. 


The greatest increase in stability again occurred with the finest grid resolution, as seen on the 120x120 
grid in figure 16. The 1-2-1 scheme reached a CFL« 20 with 12 overlaps and the 2-4-2 scheme reached a 
CFL« 18 with 12 overlaps. The 2-4-2 scheme was unstable with less than 11 overlaps, as was the 1-2-1 
scheme with less than 7 overlaps. 




(a) SBP 1-2-1 Scheme 


(b) SBP 2-4-2 Scheme 


Figure 16. CFL vs. Time step with an error tolerance of e = 5E-2 on the 120x120 grid. 


Table 3 shows the factor of increase, which is the amount of increase the OP IMEX scheme achieved as 
compared to the fully explicit scheme, for both schemes (SBP 1-2-1 and SBP 2-4-2) on the four different 
grid sizes ran with the largest controller error tolerance, e = 5E-2. The unstable results are marked with a 
The trends explained above are seen in the tables below. The tables are provided for easy comparison 
with the results from Burger’s equation. The greatest factor of increase was « 10 for the 1-2-1 scheme and 
ss 12 for the 2-4-2 scheme, both solved with the 120x120 grid. Theses results are on the same order of 
magnitude as the Burger’s equation results with the lowest viscosity, which is reasonable because the 2D 
advection equation is purely convective (no dissipation). 
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Table 3. Stability comparison of OP IMEX and explicit schemes solved on a grid with 30x30 points. 


1^=30, N y = 30 



SBP 1-2-1 

SBP 2-4-2 


CFL 

Factor of Increase 

CFL 

Factor of Increase 

Fully Explicit 

2.0 

1.0 

1.5 

1.0 

IMEX, overlap = 0 

- 

- 

- 

- 

IMEX, overlap = 1 

7.5 

3.8 

- 

- 

IMEX, overlap = 2 

7.4 

3.7 

- 

- 

IMEX, overlap = 3 

7.7 

3.9 

6.6 

4.4 

IMEX, overlap = 4 

8.0 

4.0 

7.0 

4.7 


N x =6 0, N y =6 0 



SBP 1-2-1 

SBP 2-4-2 


CFL 

Factor of Increase 

CFL 

Factor of Increase 

Fully Explicit 


2.0 

1.0 

1.5 

1.0 

IMEX, overlap 

= 3 

- 

- 

- 

- 

IMEX, overlap 

= 4 

11.6 

5.8 

- 

- 

IMEX, overlap 

= 5 

12.0 

6.0 

- 

- 

IMEX, overlap 

= 6 

12.2 

6.1 

10.6 

7.1 

IMEX, overlap 

= 7 

- 

- 

11.0 

7.3 

IMEX, overlap 

= 8 

12.6 

6.3 

11.2 

7.5 

IMEX, overlap 

= 9 

12.7 

6.4 

11.2 

7.5 


jVa;=90, Ny = 90 



SBP 1-2-1 

SBP 2-4-2 


CFL 

Factor of Increase 

CFL 

Factor of Increase 

Fully Explicit 

2.0 

1.0 

1.5 

1.0 

IMEX, overlap = 6 

- 

- 

- 

- 

IMEX, overlap = 7 

15.7 

7.9 

- 

- 

IMEX, overlap = 8 

15.8 

7.9 

- 

- 

IMEX, overlap = 9 

16.0 

8.0 

14.3 

9.5 

IMEX, overlap = 10 

16.1 

8.1 

14.4 

9.6 


JVa; = 120, Ny = 120 



SBP 1-2-1 

SBP 2-4-2 


CFL 

Factor of Increase 

CFL 

Factor of Increase 

Fully Explicit 


2.0 

1.0 

1.5 

1.0 

IMEX, overlap = 

6 

- 

- 

- 

- 

IMEX, overlap = 

7 

18.5 

9.3 

- 

- 

IMEX, overlap = 

8 

- 

- 

- 

- 

IMEX, overlap = 

9 

19.0 

9.5 

- 

- 

IMEX, overlap = 

10 

19.2 

9.6 

- 

- 

IMEX, overlap = 

11 

19.5 

9.8 

17.2 

11.5 

IMEX, overlap = 

12 

19.5 

9.8 

17.7 

11.8 
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As a sanity check on the solutions obtained with the largest time step controller tolerance, the L 2 error 
norms are presented in table 4. They are computed using the approximate solution and the exact analytic 
solution. The L 2 error norm of the explicit solution ran with the highest stable controller tolerance are also 
presented for comparison. 

The majority of the solutions have an error with the same order-of-magnitude as the explicit solution, or 
one order-of-magnitude larger than the explicit solution. There appears to be no trends with the number of 
overlaps or SBP scheme and the amount of error. The time step controller shows oscillatory behavior when 
near the stability envelop which could cause a scatter in errors reported here. When the solution ends, if it 
is at the upper end of the oscillation the error would be larger, whereas the opposite would be observed if 
it’s trending downward or decreasing the time step when the solution ends. Also, the L 2 error norms are 
expected to be larger when at the edge of stability as in these cases. This confirms that the results presented 
in this section are reliable. 


Table 4. Comparison of L 2 error norms of fully explicit and OP IMEX schemes on 120x120 grid. 


Z/2 Error Norms 


N x = 120, N v = 120 

SBP 1-2-1 

SBP 2-4-2 

€ 

Fully Explicit 

0.01081 

0.00931 

0.002 

IMEX, overlap = 7 

0.20030 

- 

0.05 

IMEX, overlap = 8 

- 

- 

0.05 

IMEX, overlap = 9 

0.03186 

- 

0.05 

IMEX, overlap = 10 

0.08419 

- 

0.05 

IMEX, overlap =11 

0.13228 

0.04079 

0.05 

IMEX, overlap = 12 

0.16849 

0.13382 

0.05 

L; 

; Error Norms 


Aj,=90, A a =90 

SBP 1-2-1 

SBP 2-4-2 

€ 

Fully Explicit 

0.03661 

0.02579 

0.0024 

IMEX, overlap = 7 

0.31106 

- 

0.05 

IMEX, overlap = 8 

0.03421 

- 

0.05 

IMEX, overlap = 9 

0.08305 

0.05441 

0.05 

IMEX, overlap = 10 

0.11561 

0.11999 

0.05 

A 

; Error Norms 


A x =60, Ay=60 

SBP 1-2-1 

SBP 2-4-2 

€ 

Fully Explicit 

0.01451 

0.05119 

0.0058 

IMEX, overlap = 4 

0.34242 

- 

0.05 

IMEX, overlap = 5 

0.08720 

- 

0.05 

IMEX, overlap = 6 

0.17650 

0.17580 

0.05 

IMEX, overlap = 7 

- 

0.30452 

0.05 

IMEX, overlap = 8 

0.33656 

0.38389 

0.05 

IMEX, overlap = 9 

0.38799 

0.06696 

0.05 


Z/2 Error Norms 


A r =30, Ay =30 

SBP 1-2-1 

SBP 2-4-2 

€ 

Fully Explicit 

0.13594 

0.03471 

0.0157 

IMEX, overlap = 0 

- 

- 

0.05 

IMEX, overlap = 1 

0.16943 

- 

0.05 

IMEX, overlap = 2 

0.14951 

- 

0.05 

IMEX, overlap = 3 

0.34547 

0.09775 

0.05 

IMEX, overlap = 4 

0.47630 

0.25975 

0.05 
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V. Design Order 


Verifying the design order is an important validation of the mathematical correctness of the 2D OP IMEX 
method. The spatial and temporal convergence tests are presented in Figure 17, showing that the OP IMEX 
methodology and the implementation are correct for the 2D advection equation. All plots are in log-log scale 
in order to more readily see the order of convergence (as a slope rather than a power). Tables containing 
the exact design order error values are shown in appendix A. 
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N N 

(a) Spatial design order of two: No overlaps (b) Spatial design order of two: Two overlaps 



N 



N 


(c) Spatial design order of three: No overlaps 


(d) Spatial design order of three: Two over- 
laps 



(e) Temporal design order of four: No over- 
laps 


0 r 



(f) Temporal design order of four: Two over- 
laps 


Figure 17. L 2 and Loo design order for SBP 1-2-1 spatial scheme with IMEX temporal advancement. 
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VI. Conclusion 


A high-order implicit-explicit multi-block time stepping method (OP IMEX) for hyperbolic PDEs was 
formulated and tested on the ID nonlinear viscous Burger’s equation and the 2D linear advection equation. 
The first and second derivatives were approximated using both second- and fourth-order accurate SBP 
operators (FD based), SBP 1-2-1 and SBP 2-4-2, respectively. 16 The semi-discrete system of ODEs is 
integrated in time using additive 6-stage Runge-Kutta schemes; 15 an explicit scheme (ERK) and an implicit 
scheme (ESDIRK), which are both fourth-order accurate in time. 

The computational domain, ID or 2D, was split into partitions or blocks. The points inside the blocks 
were implicitly solved independently of the rest of the domain (at each RK stage), while the terms coupling 
the blocks together were solved explicitly. The numerical stability of the method was evaluated by comparing 
the maximum CFL and DFL attainable while maintaining a stable solution against those of the fully explicit 
scheme. Solving the partitioned domain using the IMEX scheme alone did not increase the stability. However, 
overlapping portions of the blocks and solving the overlapped portions twice significantly increased the 
numerical stability. 

The IMEX scheme with overlapping partitions significantly reduces the CFL and DFL limitations for 
numerical stability for both the second-order (SBP 1-2-1) and the fourth-order (SBP 2-4-2) spatial operators. 
The grid resolution and viscosity studies show that the greatest increase in the stability limitations occurs 
with more viscous problems. As the viscosity is increased stability limitations decrease. As the number of 
equispaced nodes is increased, the stability limitations decrease as well. 

A. Implicit Solutions Using Direct Solving Methods 

Implicit solutions require solving the entire domain simultaneously. The large dimensionality of full 3D CFD 
simulations, which can be on the order of IE 8 - 1E9, force the use of iterative methods simply because direct 
solving methods for systems of that size are not feasible. However, the OP IMEX method presented in this 
work may provide a means for using direct solving methods in full 3D CFD. 

Nested Dissection is the optimal ordering for a direct solving method on a 3D domain, and scales as N 6 
operations for a nearest neighbor connectivity pattern. Assuming a processor can do 1E9 FLOPS (FLoating- 
point Operations Per Second) and that the 3D domain is partitioned into blocks of the size NxNxN, then 
if N « 30 the block can be inverted in ss 1 second by one CPU. If each partition is solved by an individual 
processor in a parallel fashion, direct solving methods for implicit methods used in full 3D CFD become 
feasible. 

The blocks are overlapped enough to gain the added increase in numerical stability, yet they are kept 
to a manageable size for inverting the linear matrices using direct solving methods. This newly developed 
OP IMEX method requires further evaluation, but it appears to have significant new advantages for CFD 
simulations. 


Appendix 

The ESDIRK and ERK schemes used in this work are presented below. 

Runge-Kutta Schemes 
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Fourth-Order Explicit Runge-Kutta (ERK) Scheme 


0 

0 

0 

0 

0 

0 

0 

1 

2 

1 

2 

0 

0 

0 

0 

0 

83 

250 

13861 

62500 

6889 

62500 

0 

0 

0 

0 

31 

116923316275 

2731218467317 

9408046702089 

0 

0 

0 

50 

23936840614681 

15368042101831 

11113171139209 

17 

451086348788 

2682348792572 

12662868775082 

3355817975965 

0 

0 

20 

29024286899091 

7519795681897 

11960479115383 

11060851509271 

1 

647845179188 

73281519250 

552539513391 

3354512671639 

4040 

0 

3216320057751 

8382639484533 

3454668386233 

8306763924573 

17871 

bi 

82889 

0 

15625 

69875 

2260 

1 

524892 

83664 

102672 

8211 

4 

bi 

4586570599 

29645900160 

0 

178811875 

945068544 

814220225 

1159782912 

3700637 

11593932 

61727 

225920 


Fourth-Order Singly Diagonally Implicit Runge-Kutta (ESDIRK) Scheme 


0 

0 

0 

0 

0 

0 

0 

1 

2 

1 

4 

1 

4 

0 

0 

0 

0 

83 

250 

8611 

62500 

1743 

31250 

1 

4 

0 

0 

0 

31 

5012029 

654441 

174375 

1 

0 

0 

50 

34652500 

2922500 

388108 

4 

17 

15267082809 

71443401 

730878875 

2285395 

1 

0 

20 

155376265600 

120774400 

902184768 

8070912 

4 

1 

82889 

0 

15625 

69875 

2260 

1 

524892 

83664 

102672 

8211 

4 

bi 

82889 

0 

15625 

69875 

2260 

1 

524892 

83664 

102672 

8211 

4 

k 

4586570599 

29645900160 

0 

178811875 

945068544 

814220225 

1159782912 

3700637 

11593932 

61727 

225920 
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Summation-by-parts Operators 

The SBP operators used in this work are presented here. They are obtained from Mattsson, 16 where the 
interested reader may find more detailed information on their formulation. 


SBP 1-2-1 Scheme 


The following SBP operators are first order accurate at the boundary and second order accurate in the 
interior. The discrete norm V and the discrete SBP operator T > i approximating ^ are given by 


V = dx 


V 


/-l 


V 1 =V~ 1 Q= — 
ax 


y 


v 


0 

-1 1 / 


The discrete SBP operator V 2 = V 1 {—M + BV) approximating and the boundary derivative 
operator BS are given by 


V 2 = 


dar 


/ 1-2 1 

\ 


/ 3 _o i 
2 ^ 2 

\ 

1 -2 1 


, BV = 7 - 

1 




dx 



1 

-2 1 



1 

V 1 

-2 1/ 


V 1 

-2 \) 


SBP 2~4~2 Scheme 

The following SBP operators are second order accurate at the boundary and fourth order accurate in the 
interior. The discrete norm V is given by 




V = dx 


48 


59 

48 


43 

48 


49 

48 


\ 




1 / 


V\ approximating 


given by 





/ 24 

f 17 

59 

34 

4 

17 

3 

34 

0 

0 

0 


1 

2 

0 

1 

2 

0 

0 

0 

0 

1 

4 

43 

59 

86 

0 

59 

86 

4 

43 

0 

0 

dx 

3 

98 

0 

59 

98 

0 

32 

49 

4 

49 

0 


0 

0 

1 

2 

0 

2 

1 


12 

3 

3 

12 


V 


•7 


The discrete SBP operator V 2 approximating ^2 is given by 


V 2 = V~ L (-M+BV) = —j 

dx 


2 

-5 

4 

-1 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

0 

4 

59 

110 

59 

4 

0 

0 

43 

43 

43 

43 

43 

1 

0 

59 

118 

64 

4 

0 

49 

49 

49 

49 

49 

0 

0 

1 

4 

5 

4 

1 

12 

3 

2 

3 

12 


•7 
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and the boundary derivative operator BV is given by, 


( X i ^ 

0 


BV= — 
dx 


V 


\ 


l 

3 



11 

6 


) 


Design Order Tables 

The L 2 and L ^ error norms and convergence rates for Burger’s equation and the two-dimensional advection 
equation are shown here. 

Burger’s Equation 


Table 5. Burger’s Equation: second-order spatial accuracy with no overlaps. 


OP IMEX SBP 1-2-1, No overlaps 


N 

1/2 Error 

L 2 Rate 

Lao Error 

Loo R&tG 

33 

4.79E-002 

- 

2.18E-001 

- 

63 

1.34E-002 

1.972369493 

7.83E-002 

1.580341848 

126 

3.27E-003 

2.032683149 

2.06E-002 

1.929890762 

258 

7.72E-004 

2.015284878 

5.12E-003 

1.941126999 

513 

1.94E-004 

2.006000201 

1.29E-003 

2.003397616 


OP IMEX SBP 1-2-1, One overlap 


N 

1/2 Error 

L 2 Rate 

Lac, Error 

Loo R&te 

33 

4.79E-002 

- 

2.18E-001 

- 

63 

1.34E-002 

1.972369493 

7.83E-002 

1.580341848 

126 

3.27E-003 

2.032683149 

2.06E-002 

1.929890762 

258 

7.72E-004 

2.015284878 

5.12E-003 

1.941126999 

513 

1.94E-004 

2.006000202 

1.29E-003 

2.003397617 


OP IMEX SBP 2-4-2, No overlaps 


N 

1/2 Error 

L 2 Rate 

Lao Error 

Loo R&tG 

33 

2.01E-002 

- 

8.00E-002 

- 

63 

2.45E-003 

3.255233370 

1.38E-002 

2.713134419 

126 

1.84E-004 

3.735394185 

1.22E-003 

3.502319399 

258 

1.10E-005 

3.924161809 

8.07E-005 

3.791414429 

513 

7.14E-007 

3.984038129 

5.23E-006 

3.982542337 


OP IMEX SBP 2-4-2, One overlap 


N 

L 2 Error 

L 2 Rate 

Lao Error 

Lqo R&tc 

33 

2.01E-002 

- 

8.00E-002 

- 

63 

2.05E-003 

3.530096975 

9.26E-003 

3.334727521 

126 

1.84E-004 

3.478977829 

1.22E-003 

2.922444018 

258 

1.10E-005 

3.924161797 

8.07E-005 

3.79141443 

513 

7.14E-007 

3.984033947 

5.23E-006 

3.982542623 
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Table 6. Burger’s Equation: 


OP IMEX, No overlaps 


dt 

L 2 Error 

L 2 Rate 

0.18 

5.50E-006 

- 

0.09 

3.76E-007 

3.871331123 

0.045 

2.42E-008 

3.954481746 

0.0225 

1.59E-009 

3.933310419 

0.01125 

1.37E-010 

3.535462943 

0.005625 

1.84E-011 

2.892710861 


fourth-order temporal accuracy. 


OP IMEX, One overlap 


dt 

L 2 Error 

L 2 Rate 

0.18 

5.48E-006 

- 

0.09 

3.75E-007 

3.869031493 

0.045 

2.48E-008 

3.916422375 

0.0225 

2.07E-009 

3.584320001 

0.01125 

2.80E-010 

2.88875141 

0.005625 

3.90E-011 

2.841845516 
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Two-dimensional Advection Equation 


Table 7. 2D Advection Equation: second-order spatial accuracy with no overlaps. 


OP IMEX SBP 1-2-1, No overlaps 


N X =Ny 

Li Error 

CD 

Lqo Error 

Lqo Rate 

15 

2.17E-003 

- 

5.88E-003 

- 

30 

4.94E-004 

2.134640187 

1.54E-003 

1.930163796 

60 

1.13E-004 

2.127598038 

3.60E-004 

2.098610437 

120 

2.72E-005 

2.054168695 

8.92E-005 

2.01334724 

240 

6.68E-006 

2.028252476 

2.22E-005 

2.007627234 


OP IMEX SBP 1-2-1, One overlap 


II 

Sri 

Li Error 

Li Rate 

Loo Error 

Lqo Rate 

15 

2.17E-003 

- 

5.88E-003 

- 

30 

4.94E-004 

2.134641331 

1.54E-003 

1.928978607 

60 

1.13E-004 

2.127807799 

3.60E-004 

2.100561594 

120 

2.72E-005 

2.054103722 

8.92E-005 

2.013569205 


OP IMEX SBP 2-4-2, No overlaps 


N X =Ny 

Li Error 

Li Rate 

Loo Error 

Lqo Rate 

21 

3.12E-004 

- 

1.08E-003 

- 

42 

3.77E-005 

3.051580362 

1.53E-004 

2.821834935 

84 

3.85E-006 

3.289190695 

1.98E-005 

2.950534388 

168 

4.49E-007 

3.100397842 

2.31E-006 

3.10250198 


OP IMEX SBP 2-4-2, One overlap 


II 

s? 

Li Error 

Li Rate 

Loo Error 

Lqo Rate 

21 

3.12E-004 

- 

1.08E-003 

- 

42 

3.77E-005 

3.05156826 

1.53E-004 

2.821815868 

84 

3.85E-006 

3.289191155 

1.98E-005 

2.95049937 

168 

4.49E-007 

3.100216603 

2.31E-006 

3.103500644 


Table 8. 2D Advection Equation: fourth-order temporal accuracy. 


OP IMEX, No overlaps 


dt 

L 2 Error 

L 2 Rate 

0.200000 

6.63E-003 

- 

0.100000 

8.17E-004 

3.021287254 

0.050000 

8.07E-005 

3.340077632 

0.025000 

6.40E-006 

3.656917852 

0.012500 

4.55E-007 

3.812723838 

0.006250 

3.05E-008 

3.898016031 

0.003125 

1.98E-009 

3.947159646 


OP IMEX, Two overlaps 


dt 

L 2 Error 

L 2 Rate 

0.200000 

3.36E-003 

- 

0.100000 

4.90E-004 

2.776292772 

0.050000 

5.42E-005 

3.177890808 

0.025000 

4.59E-006 

3.560298861 

0.012500 

3.35E-007 

3.775275458 

0.006250 

2.26E-008 

3.888408113 

0.003125 

1.47E-009 

3.945146214 
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